//
// Created by tong on 9/27/24.
//

#pragma once
#include <string>
#include <vector>
#include <fstream>

#include "Primitive.h"
#include "GRDECLRaw.h"

/**
 * 角点网格几何处理与角点网格转换类.
 */

class GRDECLGrid {
public:
    // 枚举类型:角点网格坐标系的类型
    enum class CoordinateSystemType { Inconclusive, RightHanded, LeftHanded };
    // 枚举类型：网格面的类型
    enum class FaceType { IFace, JFace, KFace,};
    // 枚举类型：单元中心的类型:质心、几何中心
    enum class CellCenterType { Barycenter, GeometricCenter,};
    // 定义面为顶点号的索引数组
    using Face = IntList;
    // 定义面列表为顶点号的索引的二维数组
    using FaceList = IntListList;
    // 定义面类型列表
    using FaceTagList = vector<FaceType>;


    // 顶点距离的容差
    static double tol_; //Absolute tolerance of node-coincidence
    // 是否启用性能日志输出
    static bool perfFlag_;
    // 是否Pinch collapseed网格
    static bool pinchActive_;
    // 是否输出调试文件
    static bool writeDebugFiles_;
    // 网格中心的类型
    static CellCenterType cellCenterType_;
    // 构造函数
    explicit GRDECLGrid(const GRDECLRawData& rawGRDECL);


//private:
    // GRDECL原始数据
    const GRDECLRawData& rawGRDECL_;
    // GRDECL原始网格数量(不考虑是否active)
    const int nRawCell_;
    // z方向顶点变化符号 zcornSign_ = 1 为深度严格递减
    const int zcornSign_;
    // 网格坐标系的类型
    CoordinateSystemType CoordSysType_;
    // 实际网格数量(减掉non-active网格)
    int nCell_;

    // GRDECL grid index data

    // 实际active网格到原始GRDECL网格的映射
    IntList localCellGlobalIndex_;
    //原始GRDECL网格到实际active网格的映射
    IntList globalCellLocalIndex_;
    // 每个角点位置的顶点数据
    CornerListListList cornerMap_;
    // 每个角点位置的顶点编号
    IntListListList posPointIndex_;
    // 断层错位后产生的新交点（node）数据
    IntersectIndexList intersectIndices_;
    
    /***************GRDECL网格数据***************/
    // GDECL网格的顶点坐标
    PointList nodeCoord_;
    // GDECL网格的网格面
    FaceList faces_;
    // GDECL网格的网格面类型
    FaceTagList faceTags_;
    // GDECL网格的网格面邻接关系
    IntList faceNeighbors_;
    // GDECL网格断层面索引
    IntList faultFaces_;
    // BoundingBox
    pair<Point, Point> boundingBox_;

    /***************多边形网格数据***************/
    //网格面中心
    PointList faceCenters_;
    //网格面法向量
    VecList faceAreaNormals_;
    //网格体心
    PointList cellCenters_;
    //网格体体积
    DoubleList cellVolumes_;
    //网格面索引
    IntListList cellFaces_;
    //网格顶点索引(所有顶点)
    IntListList cellNodes_;
    //网格顶点索引(先K顶面后K底面)
    IntListList cellKNodes_;
    //IntListList cellEdgeFaces_;
    Map<int,int> faceKValue_;

    /***************多边形网格属性数据***************/
    //杨氏模量
    DoubleList YoungsMod_;
    //泊松比
    DoubleList Poisson_;
    //
    DoubleList Poro_;
    //
    DoubleList RockDens_;

    /***************四面体分解后网格数据***************/
    //四面体网格顶点坐标
    PointList allTetCoord_;
    //四面体网格单元顶点索引
    IntListList allTetCellNodes_;
    //四面体网格单元数据
    int tetReserviorCellNum_;
    IntList tetPointData_;
    IntList tetCellData_;
    PointList tetCellCenter_;
    DoubleList tetCellVolume_;
    //
    DoubleList tetYoungsMod_;
    DoubleList tetPoisson_;
    DoubleList tetPoro_;
    DoubleList tetRockDens_;

    /***************六面体分解后网格数据***************/
    //六面体网格顶点坐标
    PointList allHexCoord_;
    //六面体网格单元顶点索引
    IntListList allHexCellNodes_;
    //六面体网格单元数据
    int hexReserviorCellNum_;
    IntList hexPointData_;
    IntList hexCellData_;
    PointList hexCellCenter_;
    DoubleList hexCellVolume_;
    //
    DoubleList hexYoungsMod_;
    DoubleList hexPoisson_;
    DoubleList hexPoro_;
    DoubleList hexRockDens_;

    
    /***************四面体-金字塔-六面体混合分解后网格数据***************/
    //混合网格顶点坐标
    PointList decomposedCoord_;
    //混合网格单元顶点索引
    IntListList decomposedCellNodes_;
    //混合网格单元数据
    int decompReservoirCellNum_;
    IntList decomposedPointData_;
    IntList decomposedCellData_;
    PointList decomposedCellCenter_;
    DoubleList decomposedCellVolume_; 
    //
    DoubleList decomposedYoungsMod_;
    DoubleList decomposedPoisson_;
    DoubleList decomposedPoro_;
    DoubleList decomposedRockDens_;


    /***************延拓数据***************/
    // 提取的原始地层面和边界面
    PointList extendPoints_;
    IntSet IMinFace_;
    IntSet IMaxFace_;
    IntSet JMinFace_;
    IntSet JMaxFace_;
    IntSet KMinFace_;
    IntSet KMaxFace_;
    // 提取的原始断层面
    //IntList faultCoordFace_;
    // 延拓网格数据所需边界数据
    IntListList IMinFaceNodes_;
    IntListList IMaxFaceNodes_;
    IntListList JMinFaceNodes_;
    IntListList JMaxFaceNodes_;
    IntListList KMinFaceNodes_;
    IntListList KMaxFaceNodes_;





    /*************************几何处理成员函数*************************/
    //获取角点位置唯一顶点坐标
    void getNodeCoord();
    //获取垂直网格面
    void getVerticalFace();
    //获取水平网格面
    void getHorizontalFaces();
    //计算多边形网格数据
    void getNonSlipFaultFace();
    //计算多边形网格断层数据
    void getCommNormFaultFace();
    //计算新的相交顶点
    void calcIntersectPoints();
    //重排GRDECL数据
    void reorderGRDECLData();
    //计算多边形网格数据
    void calcPolyGridData();
    //全四面体网格分解
    void decomposeAllTetCell();
    void decomposeAllTetExtendCell();
    //全六面体网格分解
    void decomposeAllHexCell();
    //四面体-金字塔-六面体混合网格分解
    void decomposeCell();
    // 生成地层和延拓界面
    void genExtendCell(PointList &nodeCoord,
                       IntListList &cellNodes,
                       IntList &pointData,
                       IntList &cellData);
    // 生成地层和延拓界面
    void genHybridExtendCell();
    void decomposeHybridExtendCell();
    void genHexExtendCell();
    void genTetExtendCell();

    static double extendZMax_;
    static double extendZMin_;
    static double extendXFactor_;
    static double extendYFactor_;
    static int XExtendNum_;
    static int YExtendNum_;
    static int ZExtendNum_;

    // 在Pillar上获取内部垂直网格面
    // 输入：两个Pillar对上角点位置的索引
    // 输出: 每个垂直网格面两侧单元的K方向索引
    IntList findVerticalFace(
        const IntList& a1,
        const IntList& a2,
        const IntList& b1, 
        const IntList& b2);

    // 获取网格坐标系统类型
    // 输出:坐标系类型（左手坐标系/右手坐标系/无法判断）
    CoordinateSystemType getCoordSysType() const;
    
    // 获取ZCORN数据符号（1为递增、-1为递减）
    // 输出:ZCORN数据递增或递减符号
    int getZCORNSign() const; 

    // 获取网格垂向尺寸
    // 输入：网格单元线性索引，网格八个角点位置序号索引
    // 输出: 网格单元垂向尺寸
    inline double vertSize(const size_t c, const size_t off[8]) const noexcept;

    // 判断是否为垂直方向上相邻
    // 输入：两个网格单元线性索引
    // 输出: 是否为垂直方向上相邻
    inline bool isVerticalNeighbor(const int c1, const int c2) const noexcept;

    // 获取网格单元的线性表索引
    // 输入：网格单元IJK索引
    // 输出: 网格单元线性表索引
    inline int cellLinearIndex(const int I, const int J,const int K) const noexcept;
    // // 判断是否为K方向上相邻
    // inline bool isKNeighbor(const int c1,const int c2) const noexcept;

    // 输出VTK格式网格几何数据
    // 输入：输出vtk网格名
    void vtkCellWriter(const std::string &outputPath) const;

    // 输出VTK格式网格体数据
    // 输入：输出vtk网格名
    void vtkCellDataWriter(const std::string &outputPath,const DoubleList& data,const DoubleList& data2) const;

    // 输出全六面体分解后网格几何
    // 输入：输出vtk网格名
    void vtkHexCellWriter(const std::string &outputPath);

    // 输出全四面体分解后网格几何
    // 输入：输出vtk网格名
    void vtkTetCellWriter(const std::string &outputPath);

    // 输出混合分解后网格几何
    // 输入：输出vtk网格名
    void vtkDecomposeCellWriter(const std::string &outputPath) const;
    
    // 输出断层面几何数据
    // 输入：输出vtk网格名
    void vtkFaultFaceWriter(const std::string &outputPath) const;

    // 输出断层面几何数据
    // 输入：输出vtk网格名
    void vtkBoundaryWriter(const std::string &outputPath) const;

    void vtkKMaxWriter(const std::string &outputPath) const;

    // 输出断层面几何数据
    // 输入：输出vtk网格名
    void vtkExtendedCell(const std::string &outputPath);

    // 输出断层面几何数据
    // 输入：输出vtk网格名
    static void vtkExtendedCell(
        const std::string &outputPath,
        const PointList &nodeCoord,
        const IntListList &cellNodes,
        const IntList &pointData,
        const IntList &cellData);

    // 输出断层面几何数据
    // 输入：输出vtk网格名
    static void vtkExtendedCellData(
        const std::string &outputPath,
        const PointList& nodeCoord,
        const IntListList& cellNodes,
        const IntList& pointData,
        const IntList& cellData,
        const DoubleList& youngsMod,
        const DoubleList& poisson,
        const DoubleList& poro,
        const DoubleList& rockDens);

    // 输出赛题要求的属性数据
    // 输入：输出文件名
    void writeRequiredData(const std::string& outputPath);

    void scaleCellData(DoubleList& data, const double scale) const;

    static void assignCellData(
        const IntList& cellTypes,
        DoubleList& data,
        const double overValue,
        const double underValue,
        const double sideValue,
        const double overSideValue,
        const double underSideValue);

    void calcExtendCellCenterAndVolume(
        const PointList& nodeCoord,
        const IntListList& cellNodes,
        PointList& cellCenter,
        DoubleList& cellVolume) const;

    // 计算守恒值
    static double calcConservation(const DoubleList& value,const DoubleList& coeff,const size_t dataSize);
    // 缩放数据
    static void scaleData(DoubleList& data,const size_t dataSize, const double scale);
    // 提取地层面
    static void extractLayer(IntList& pointData,const IntList& layers);

    // 释放部分冗余内存 for 赛题1
    void freeGrdeclData();
    void freeHexData();
    void freeTetData();
    void freeDecomposedData();
};



// inline function for GRDECLGrid
double GRDECLGrid::vertSize(const size_t c, const size_t off[8]) const noexcept
{
    size_t i = c % rawGRDECL_.dimI;
    size_t cd = c / rawGRDECL_.dimI;
    size_t j = cd % rawGRDECL_.dimJ;
    size_t k = cd / rawGRDECL_.dimJ;
    const double *zcorn = rawGRDECL_.zcoord.data();
    size_t start = (2 * i) + (2 * rawGRDECL_.dimI) * ((2 * j) + (2 * rawGRDECL_.dimJ) * (2 * k));
    double dz = 0.0;
    for (k = 0; (!(fabs(dz) > 0)) && (k < 4); k++)
    {
        dz = zcorn[start + off[k + 4]] - zcorn[start + off[k]];
    }
    return dz;
}

bool GRDECLGrid::isVerticalNeighbor(const int c1, const int c2) const noexcept
{
    int k1 = c1 / rawGRDECL_.dimI / rawGRDECL_.dimJ;
    int k2 = c2 / rawGRDECL_.dimI / rawGRDECL_.dimJ;
    return std::abs(k2 - k1) == 1;
}

int GRDECLGrid::cellLinearIndex(const int I, const int J, const int K) const noexcept
{
    return I + rawGRDECL_.dimI * (J + rawGRDECL_.dimJ * K);
}

// bool GRDECLGrid::isKNeighbor(const int c1, const int c2) const noexcept
// {
//     const int k1 = c1 / rawGRDECL_.dimI / rawGRDECL_.dimJ;
//     const int k2 = c2 / rawGRDECL_.dimI / rawGRDECL_.dimJ;
//     return std::abs(k1 - k2) == 1;
// }
